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Abstract 

In this work, transverse shear warping functions for an equivalent single layer plate model are formulated from 
a variational approach. The part of the strain energy which involves the shear phenomenon is expressed in 
function of the warping functions and their derivatives. The variational calculus leads to a differential system 
of equations which warping functions must verify. Solving this system requires the choice of values for the 
(global) shear strains and their derivatives. A particular choice, which is justified for cross-ply laminates, leads 
to excellent results. 

For single layer isotropic and orthotropic plates, an analytical expression of the warping functions is given. 
They involve hyperbolic trigonometric functions. They differ from the z — A/3z 3 Reddy's formula which has 
been found to be a limit of present warping functions for isotropic and moderately thick plates. When the h/L 
and/or the G13/E1 ratios significantly differ from those of isotropic and moderately thick plates, a difference 
between present warping functions and Reddy's formula can be observed, even for the isotropic single layer 
plate. Finite element simulations agree perfectly with the present warping functions in these cases. 

For multilayer cross-ply configurations, the warping functions are determined using a semi-analytical proce- 
dure. They have been compared to results of 3D finite element simulations. They are in excellent agreement. 

For angle-ply laminates, the above process gives warping functions that seem to have relevant shapes, even 
if the choice for global shear values cannot be justified in this case. No finite element comparison has been 
presented at this time because it is difficult to propose boundary conditions and prescribed load that permit to 
isolate the shear phenomenon. 

Keywords: Multilayered, anisotropic, plate, laminate, warping functions, transverse shear 



1. Introduction 

The need to take into account the transverse shear behavior in plate theories has emerged when authors 
began to focus on thick homogeneous plates or on inhomogeneous plates. For homogeneous plates, the works of 
Reissner [1], Uflyand [2] and Mindlin [3] are considered as historical references. Industrial composite materials 
have appeared in the end of the 1930s. They can easily be used into multilayered structures, with some benefits. 
For these particular anisotropic materials, and also for multilayered plates, even for those constituted of isotropic 
plies, it has been rapidly evident that the weak and/or inhomogeneous transverse shear stiffness permit the shear 
deformation to significantly appear. Authors have then proposed models that fit this behavior. Among them, 
we can cite the early works of Lekhnitskii [4] and Ambartsumyan [5]. However, the one-layered models based 
on the Mindlin-Reissner model had been adapted to multilayered structures, giving the nowadays called first 
order shear deformation theory (FOSDT). This theory is particularly efficient when shear correction factors are 
used. Those shear correction factors take into account the variations of the transverse shear stresses through 
the thickness. In what has become a reference work, Whitney [6] gave a method for their computation. Other 
authors have later proposed different ways to compute them [7, 8] . 

Since these works, static and dynamic behavior of anisotropic multilayered plates have been extensively 
studied, as it can be seen in Carrera's and Khandan's review articles [9, 10]. Two important classes of plate 
theories that appear in such reviews are commonly called: the equivalent single layer theories (ESL) and the 
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layerwise theories (LT). The latter contains theories that are based on hypothesis which are made in each layer, 
and lead to a number of variables which -generally- increases with the number of layers. Some of them however 
propose an elimination process which lead to a number of variables which does not depend on the number of 
layers. On the other hand, the first class contains theories that are based on variables defined in a reference 
plane, and generally propose polynomial developments of fields over z. However, some of them, the higher 
order theories propose advanced kinematic variations through the thickness, for example introducing warping 
functions which generally are built from the ply sequence. Accounting to these recent developments, it seems 
that there is a very thin line between these two class. 

For all classes combined, the transverse shear strain variations through the thickness have recently received 
many attention, as it can be seen in references [11-15]. For the ESL class, the more recent and pertinent 
theories have been built using warping functions [11, 16]. These warping functions introduce the through-the- 
thickness (z) transverse shear strain variations directly in the layer displacement field. In these models, the first 
in-plane (x and y) spatial derivatives of the (global) transverse shear strains are coupled with the membrane 
and bending terms, leading to a more complex stiffness matrix of size 10 x 10 compared to the classical 6x6 
stiffness matrix of the FOSDT model. The behavior of these models are mainly determined by the choice of 
the warping functions. This is an interesting point because warping functions can be changed and/or adapted 
on the fly to fit particular requirement like special materials (adjacent plies with high Young's modulus ratio, 
viscoelastic behavior, functionally graded materials...) and/or in-plane inhomogeneous structures (damped 
structures with patches, structures with inserts, windows...). Indeed, the change of warping functions only 
requires a new computation of the stiffness and mass matrices of the model, which consists on evaluations of 
integrals over z. The process can therefore be completely automated. All the variables of the model remains 
the same, and keep the traditional meaning of ESL variables: in-plane displacements, transverse displacement, 
rotations and transverse shear, all these quantities being associated with a given reference plane. 

There are several techniques that can be used to formulate warping functions: 

1. One can try to built them by integration of the third equilibrium equation. Functions, which can be 
assimilated to warping functions, appear during the shear correction factor computation [6]. 

2. One can also propose to search them into a reduced function family, like piecewise constant or quadratic 
functions of z, and then try to find the unknowns coefficients by means of kinematic and static conditions 
at interfaces [8]. An early work of Sun & Whitney [17], which has been later developed and enhanced 
to more general problems [18], has been entirely reformulated in reference [16] with the help of piecewise 
constant warping functions. 

3. They also can be built from an energetic variational approach [11]. In these approaches, it is also possible 
to restrain the functional space to constant or quadratic functions. 

In the present study, warping functions are obtained from an energetic variational approach. There is no 
a priori hypothesis on their shapes. The main equations of the associated plate model [16] are recalled here. 
Next, the part of the strain energy that involves the warping functions is isolated, and reformulated in order to 
make the variational calculus. It leads to a differential system of equations that warping functions must verify. 
The global values of the transverse shear strains and their spatial derivatives must be fixed in order to make 
the system solvable. A particular choice, which is justified for cross-ply laminates, gives excellent results. For 
angle-ply laminates, this choice is kept, but further work must be done to precise if this choice is relevant or if 
it exists a better way to do it. Then, an analytical solution is given, involving unknown constants for each layer 
of the multilayered structure. For a one-layer orthotropic plate, the constants are easy to determine and then 
the solution can be given in an analytical form. For multilayered plates, the solution requires the determination 
of these unknown constants which are related to continuity relations prescribed at the interfaces. 

2. Associated plate model 

In this section, the main equations of the associated plate model of reference [16] is presented. In the 
following, Greek subscripts take values 1 or 2 and Latin subscripts take values 1, 2 or 3. The Einstein's 
summation convention is used for subscripts only. The comma used as a subscript index means the partial 
derivative with respect to the directions corresponding to the following indexes. 
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Figure 1: Geometrical parameters of the undeformed laminate on the left side, and its deformed shape with corresponding quantities 
on the right side. 

2.1. Laminate definition 

The laminate is composed of n layers located between —h/2 and h/2, where h is the total height. The 
reference plane is taken as the z = plane. Figure 1 helps to visualize the following definitions: 

- z l is the elevation/offsetting of the middle plane of the layer I 

- the £th layer is located between elevations £ f_1 and Cf 

With these definitions: 

- there are n parameters z l , 

- there are n + 1 parameters £* with i taking values from to n, 

- the thickness of the layer £ is h = ( e — C 1 ^ 1 

2.2. Displacement field 

The kinematic assumptions of the present formulation are: 

f u a (x,y,z) = u° a (x,y)~ zw° a (x,y)+ip a p(z)j° p3 (x,y) _ 
1 u 3 (x,y,z) = w°(x,y) 

where u° a (x,y), 7^3(2;, y), and w°(x,y) are respectively the in-plane displacements, the engineering transverse 
shear strains 1 , and the transverse displacement at the reference plane, and ip a p{z) are four warping functions 
that will be determined. 



^^Note that if the reference plane coincides with an interface, this definition is not correct, but we shall see that it is always 
possible to define equivalent transverse strains. 
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2.3. Strain field 

Let us compute the strain field from equations (1): 



e a8 (x,y,z) = e° aB (x,y) - zw° aB (x : y) + - {ip a ^{z)^ B {x,y) + (p Bj (z)tf 3}Ct (x, y)) 
e a3 (x,y,z) = ^p' a p{z)^{x,y) 
£33(2, y,z) = 



(2a) 
(2b) 
(2c) 



Note that e° aBl w aj3 , 7^3 B and 7^3 form a set of 12 independent generalized strains to consider in this model. 



2.4. Stress field 

The strains obtained in formula (2) permit to compute the stresses: 

a aB (x, y, z) = Q a/37 s(z) (e° s (x, y) - zw° 7 g(x, y) + ip lf ,(z)^ l3S {x, y)) 
a a3 (x, y, z) = C a3 p 3 ip' Btl {z)il 3 (x, y) 
0-33(3;, V,z) = 



(3a) 
(3b) 
(3c) 



where Q aBl s are the reduced generalized plane strain stiffnesses. The vanishing of the term | of equation (2a) 
into equation (3a) is not evident. Let us demonstrate it, omitting the x, y, and z coordinates: 



2<9a/37<5 (V7m7°3,« + ^7°3, 7 ) = ^ <9"07<5 W7°3,5 + -jQafi^S^l^^ 

= 2 < 3«/97'5V :, 7At7°3,<5 + 2 Qa^ 7 <^V7°3,<5 
= Qa/3 7 5</' 7AI 7°3, 1 5 

The vanishing of the term \ of equation (2b) into equation (3b) is due to another reason: 

Oa3 = C a3 p 3 tp 3 + C a33B e 3B = 2C a3B3 e B3 = C a3B3 (z)(p' By (z)j l3 (x, y) 



(4) 



(5) 



2.5. Kinematic and static assumptions of the model 

The model requires the continuity of in-plane displacements and transverse shear stresses at each of the n—1 
interfaces: 



lim u a (x,y,z)= lim u a (x, y, z) 
lim o- a3 (x,y, z) = lim a a3 (x,y,z) 



(6a) 
(6b) 



for £ £ [1 ... n — 1] . These conditions (6a) and (6b) can now be written in terms of conditions on the warping 
functions: 



lim <p aB (z) = lim ip aB (z) 
z->C - z->C + 

lim C a3B3 (z)cp' B (z) = lim C a3B3 {z)ip' g Jz) 

2->C z^Q t + 



(7a) 
(7b) 



Equations (7a) represent the continuity of the warping functions at the n—1 interfaces whereas equations (7b) 
represent "jump" conditions on their derivatives. 
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2.6. Strain energy 

It is possible to compute the strain energy density S J = l/2e,j<xy from formulas (2) and (3) and integrate 
it over the thickness to obtain a strain energy surface density J(x,y): 



f h/2 




/ tijVijdz 


-h/2 




rh/2 




'-h/2 




fh/2 
/ 




'-h/2 





I rh/2 

7, / (tafiVufi + 2e a3 cr Q 3 + 6330-33) dz 

1 J -h/2 



dz 



[(4/3 - zw ° a p + <Pa-y(zh° 3 ,p) a a/3 + ip' af) (z)jj 3 a a3 \ dz 
x and y have been omitted in the right hand side, for clarity. 



(8) 



2. 7. Generalized forces 

The strain energy can also be written 

J = 2 Kp n ^ ~ ™%M a p + j% P^ + *f pa Q p ] 

naturally introducing the following quantities which are the generalized forces, 

rh/2 

{N a/3 , M a0 , P 1 p} = / {l,z,y ai {z)}<j a p{z)dz 

J -h/2 
h/2 



Qp = / ip' ap (z)a , 3 (z)dz 

J —h/2 



(9) 

(10a) 
(10b) 



each associated with a corresponding generalized displacement in the strain energy formula (9). 

N a p and M a p are respectively the classical plate membrane forces and bending moments, and P a p and Q a 
are special moments associated with the warping functions, i. e. associated with the transverse shear behavior. 
Note that P a p ^ Pp a in the general case, leading to a set of 12 generalized forces. 

2.8. Generalized stiffnesses 

Let us calculate these generalized forces with the help of equations (3) and (10). They are (some explanations 
are given below): 

N a p = Aa^s^s + -B Q( 3 7< 5(— w° 7l5 ) + E al 3^s7%,s (H- a ) 
M a/3 = B aj3lS e° lS + D aPl s(-w° lS ) + Fafi^^s (lib) 

P aj3 = E lSa/3 €" s + F^Sap^W ^) + G aj3fi s7,i3,5 (11°) 
Qa = H aml % (lid) 

where the following generalized stiffnesses have been introduced: 



{^4ct/37<5j BaP^&i DapySi Paf3fiSy Fafiii&i G^fj^s} — I Qa/3js{^T Z 7 Z > l P~fV, 

(z),zip 7l i(z), tp av (z)cp lfl (z)}dz (12a) 



C 



H a3/ 33 = / ip / ia (z)C^ 3 s3f' S0 (z)dz 



(12b) 



The N a/ 3 and M a p computation is straightforward but special attention must be paid to the calculus of the 
P a /3, revealing some uncommon symmetries: 



P a /3 = [ ip m {z)a^{z)di 
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c 



(13) 



In this last expression, the EjSa/3 and Fys a p are identified with the help of the major symmetry of the Q^p^siz) 
tensor. 

The A, B and D tensors inherit the symmetries of Hooke's tensor, a symmetry for each pair of indexes, 
called the minor symmetries and the major symmetry which permits the swap of the two pairs of indexes, this 
last one being related to the existence of a strain energy. The E and F tensors loose the symmetry on the last 
pair of indexes, forcing the major symmetry to disappear. The G tensor looses the symmetry on the two pairs 
of indexes, but keeps the major symmetry: 



'for the A, B, D tensors: A 
for the E, F tensors: E 
for the G tensor: 



(14) 



G/3 a jS 7^ Ga/3jS — G.yS a /3 ^ Gs~f/3a 



Hence there are 6 independent components for A, B and D, 12 for E and F, 10 for G, and 3 for H. So this 
plate model has a total of 55 independent stiffness coefficients in the most general case. 

2.9. Static laminate behavior 

This section aims to present the static 2 laminate behavior in a matrix form. The generalized forces are set, 
by type, into vectors: 













-I 




M 




\- 




N 12j 




M l2j 






Q 



Qi 
Q2 



(15) 



and the same is done for the corresponding generalized strains: 



=11 



=12. 



-2w° 15 
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7? 3 

7 2 °3 



(16) 



Generalized forces are linked with the generalized strains by the 10 x 10 and 2x2 following behavior matrices: 

fN] [A B El [e] 

M = B D F \ K \ {Q} = [H]{ 7 } (17) 

[pj [e t f t gJ |rj 

3. Warping functions 

3.1. Variational formulation 

Replacing the generalized forces in the strain energy of formula (9) by means of formula (17), and discarding 
all the quantities which do not depend on the shear phenomenon leads to: 



1 



J s = - [r T (2E T e + 2F t k + Gr) + 7 T H 7 ] 



(18) 



2 In reference [16], the dynamic behavior is also developed, leading to generalized mass terms and motion equations. 
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We shall develop in the following lines, a matrix form of stiffnesses definitions of formula (12). Introducing the 
following matrices: 



C m — 



Qim Q1122 Q1112 

Q1122 Q2222 Qtiyi 

Q\\T2 Q2212 Ql212 



Cs = 



C2323 C2313 
C2313 C1313 



ip n ipu 

lfi 2 2 ¥>21 

ip 2 \ ipi2 ¥22 



(19) 



z being omitted for clarity, one can verify that: 



E : 



/h/2 ph/2 rh/2 rh/2 

C m *dz , F = / C m *zdz , G= / # T C m *dz , H = / * /T C s *'dz (20) 
-h/2 J-h/2 J-h/2 J-h/2 



Now, the shear contribution to strain energy can be written: 

*h/2 
I -h/2 



1 rh/2 

J s = - [T T (2* T C m e + 2z* T C m « + * T C m *r) + 7 T *' T C S *'7] dz 

1 J-h/2 



(21) 



For computational reasons, it is easier to search for continuous functions. This is the reason why we shall 
consider in the following, the "stress warping functions" 



*p a p(z) = C a3l s 3 {z)<pp 7 (z) 



(22) 



which are four functions of class C 1 . For the following calculations, we need to arrange these four functions and 
their derivative in 4 x 1 vectors: 

'ipu 



> and *' 




(23) 



Let us consider the matrices: 
, 



Z m = 4 



and: 



Si3137l3,l 
52313713,2 



5'l3237l3, 
5 , 23237l3,: 



<Sl323723,l 
"5*2323723,2 



<Sl313723,l 
£2313723,2 



.'S'l3137l3,2 + "52313713, 1 5 2 3237l3,l + 51323713,2 5i 32 3723,2 + 52323723,1 5 2 313723,1 + 5l313723,: 



It can be shown that: 



5l 3 137l3 5i3237l3 5i323723 5i3i372 3 
52313713 52323713 52323723 52313723. 



Z m * = *r and Z s *' = 



(24) 
(25) 

(26) 



Now it is possible to reformulate the shear contribution to strain energy in a manner that the warping functions 
are the variables of a quadratic form: 



J s = l [* T (2Z m C me + 2^Z m C m K + Z m C m Z m *)+*' T ZJC s Z s vI/']d2 

Z J-h/2 



I -h/2 

The variationnal calculs leads to 

rh/2 



6J S = f [5* T (Z m C m e + zZ m C m « + Z m C m Z m *) + 5*' T Zj C S Z S *'] dz 

J-h/2 



(27) 



(28) 



Note that two factor 2 had appeared during the derivation of the third and fourth terms, and then vanished 
with all others during the simplification. Integrating by parts the fourth term leads to: 



rh/2 

SJ S = / [<5* T (Z m C m e + zZ m C m K + Z m C m Z m *) - 5V T ZjC s Z s *"] dz + [<5* T Zjc s Z s *'] 

J-h/2 



,1 h/2 

h/2 



(29) 
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The last term is null if there is no prescribed shear stresses at the top and the bottom of the laminate. For 
concision, we shall denote: 



H z — ZjC s Z s ; G z — Z^C m Z m ; F z — zZ^C m ; E z — Z^C m (30) 
hence the differential system of equations that must verify the \I/ functions is written: 

H Z *" = G Z * + F z « + E z e (31) 

3.2. Solving the system 

If the matrices H z and G z have been invertible, the system (31) would have admit a general solution of the 
form 

* = C ie ^ lG ^ /2z + C 2 e- (iI *- lG ^ /2 * - G Z - X F Z « - G, 1 ^ (32) 

where Ci and C2 are two 1x4 vectors of constants defined for each layer, which must be adapted to make the 
warping functions and their derivatives verify conditions on top and bottom of the laminate, at the reference 
plane, and at the interfaces (detail is given below). Unfortunately, for the most general case, H z is of rank 2, G z 
is of rank 3 and they are both of dimensions 4x4. One can try to use generalized inverses, but the calculation 
is complicated and has not be found to give good results. The matrix H z is written: 



H z = 4 



$1313(713) S 2 313{ji 3 ) 2 52313713723 S 'l3137l3723 

$1323(713) $2323(713) $2323713723 <Sl3237l3723 

5*1323713723 $2323713723 $2323(723) $1323 (723) 

.5l3137l3723 $2313713723 $2 3 13(723) 2 5i3i 3 (72 3 ) 2 



(33) 



The general solution depends on values of 75*3 and 723, on values of their derivatives 7J3 1 , 723 2 , 713 21 723 1 
present in the matrix G z and also on values of k"ij K 22> k ?2j e ii> e 22i e i2- 

It is possible to compute locally these values, starting with a classical plate model, and then to compute the 
warping functions to build the refined model, which will give new values for these quantities, and so on, with 
no guarantee of convergence. 

Another approach consists on giving chosen values of the shear strains and their derivatives. For example, 
the choice 7 13 = 1, 7^3 = 0, 7 13 1 = —ir/L x , 723 2 = 0, 7J3 2 = 0, 723 1 = will give a reduced system of size 
2x2, only involving regular matrices, from which we can compute ^11 aid '021- The choice of —Tt/L x for the 
7i3 1 vahie will be discussed later. The corresponding choice for the y direction will give ^22 an d tpi2- 

3.3. Analytical solution for a single layer isotropic plate 

For an isotropic homogeneous plate, the system (31) reduces to two similar scalar differential equations 
whose solution is: 

■011 (z) = I a\z sinh(a 2 .z)) (34) 

a l — 1 \ 0-2 J 

ai=cosh(a 2 ^J ; a 2 = 7^ li \/^ ' ?>:,] 



where: 




The other equation gives ^22 (z) which can be obtained replacing 7J3 1 by 7 23 2 and 7° 3 by 7° 3 in the above 
formulas. 

The corresponding (strain) warping functions ip a p are obtained dividing the above function by G. In addition, 
it is possible to reduce the warping functions in a manner that they take values only between —1/2 and 1/2 
applying the rule <Pii{z) = \(pu{hz). With this choice, we obtain: 

V11 ( z ) = — ~T ( aiZ ~T siml i a 2hz) ] (36) 

a\ — 1 \ 02ft / 

The choice of the values of 7° 3 a and 7° 3 is done considering the Navier's or Levy's procedures. In these well 
known bending problems, the plate quantities (pressure, transverse displacements, transverse shears, etc.) are 
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taken as trigonometric functions. For example for a rectangular plate with < x < L x and < y < L y the 
choice lead to 713(0;, y) — 7 1 ^ ax cos(ttx/L x ) sm(ny / L y ). The derivative of this function with respect to x can 
be written 713,1(2;,?/) = 7i3 a f sm(nx/L x ) s'm(iry/L y ) that permits to write jf^f — —Tr/Lxjf^. It is the reason 
why we choose the values 7° 3 = 1, 7J3 x = —ir/L x . In addition, for an isotropic material, Q/G = 2/(1 — nu), so 
one can calculate a\ and 02 to obtain the solution which only depends on v and h/L x : 




It is interesting to compute a series expansion of this function: 

tp^z) = z + d 3 z 3 + d 5 z 5 + 0(z 6 ) (38) 

Table 1 gives values of these coefficients for v = 0.3 and h/L x G {0.0001, 0.001, 0.01, 0.1, 0.2, 0.4, 0.8}. It permits 
to evaluate the difference between this function and the Reddy's formula z — 4/3z 3 . Some of the corresponding 
warping function are also plotted and compared to finite elements simulations in section 4, figure 2. 

h/L 0.0001 0.001 0.01 0.1 0.2 0.4 0.8 

~d 3 -1.3333332 -1.3333325 -1.33326 T326 T302 T215 -0.9277 
d 5 -1.88010" 8 -1.88010" 6 -1.880KT 4 -0.01869 -0.07345 -0.2740 -0.8371 

Table 1: Coefficients and (% of the series expansion of 1/3^(2:) for the isotropic material. 



3.4- Analytical solution for a single layer orthotropic plate 

The orthotropic single layer case is interesting because, as we shall see later, it permit to set arbitrary 
small values for the transverse shear modulus. Warping functions that differ strongly from Reddy's formula can 
appear in these cases. There is another interest: the orthotropic layer can be put in an off-axis configuration, 
that lead to non null functions (pi2 and if2i- We shall see that, due to the tensorial nature of the tpa/3, it is 
possible to give exact formulas for the warping functions of an off-axis single layer orthotropic laminate. 

For an on-axis orthotropic homogeneous plate, the system (31) reduces to two uncoupled 2x2 matrix 
differential equations. In addition, the matrix H Z _1 G Z is diagonal. When the interface and top/bottom/middle 
conditions are enforced (see explanations in section 3.5), the V'2i(^) (for the first system) and the i\)n(z) (for 
the second system) are found to be null functions. Hence, the solution is: 



^ll(z) = ( aiz — sinh (a 2 z) 

01 — 1 \ a 2 

ip22(z) = ^ 23 ( a 3 z — sinh (a 4 z) 

a 3 ~l \ a 4 



(39) 
(40) 



where: 



ai = cosh 0,2 



03 = cosh I a.4 



«2 



04 



713,1 
7l°3 
723,2 



Gi 



G 



2:5 



(41) 



(42) 



As for the isotropic case, the warping functions ipn(z) and <p22(z) are obtained from ip±i(z) and ip22( z ) by 
division by G13 and G23 respectively. It is also possible to have reduced warping functions with z varying from 
— 1/2 to 1/2, in order to compare them for different ratios h/L. The choice for ratios values of 7^3 1 /^fi 3 = —n/L x 
and 7232/723 = ^^/Ly is also applied. Although it is not done in this paper, it is then possible to study 
rectangular plates. As we shall see later, three studies are done for this case, one with various ratios h/L, 
another with various ratios G\ 3 /E\ (we focus on the x direction), and an off-axis study. 
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For this last case, due to their tensor nature, warping functions tp a p must obey to classical coordinate 
transformation formulas. Hence, the warping functions of a 0-oriented single ply can be calculated from the 
warping functions of the on-axis one: 



where c = cos(0) and s — sin(0). Note that if the material have a square symmetry and if the height to length 
ratios are the same in each direction, </?5i = ¥>22 an d then the ip® 2 and y 2 i functions are null, but if it is not the 
case, they are not null. However, the generic algorithm, described in the following section, gives functions ip e 12 
and (fi2\ that are similar in shape to those of formulas (43) but do not present the symmetry which is evident in 
the above formulas. This must be studied in the future. The functions given by the code are shown in figure 5. 

3.5. Multilayered cross-ply orthotropic plates 

As it was told before, for on-axis orthotropic plies, the system (31) reduces to two uncoupled 2x2 matrix 
differential equations. Hence the system is split into two subsystems for which the solution (32) can be computed. 
The first subsystem consists of the first two lines of system (31), it permits to compute ipu and ip2i- 

There are two unknown constants per layer in Ci and two in C2- There are also three global constants in k 
and three in e. The warping functions must verify 2 conditions at the reference plane </?ai(0) = 0. The derivatives 
of the warping functions must verify 4 conditions at the top and the bottom of the plate ip' al (±h/2) = and 
2 conditions at the reference plane tp' a i{0) = <5 Ql . These conditions can be easily transformed into conditions 
for the ip stress warping functions. Hence, for a single layer plate, there are 2 + 2 + 3 + 3 = 10 constants to 
determine with 4 + 2 + 2 = 8 conditions. If the unknowns are determined using a least square solving method, 
the K12 and £12 constant are found to be null, which is coherent with the on-axis orthotropic behavior. Hence 
we have the choice to eliminate them and solve a square system, or to keep them and use a least square method. 

For multilayered plates, the warping functions ip a \ must be continuous at each interface, which gives 2 
conditions per interface, and the derivatives of the stress warping functions ip' al must also be continuous, which 
gives 2 more conditions per interface. Hence, for a n-layered plate, there are 2n + 2n + 3 + 3+ = 4n + 6 unknown 
constants and 4 + 2 + 2 + 2 (n — 1) + 2(n — 1) = 4n + 4 equations. The same remark as above, concerning the 
solving procedure, applies here. 

The second subsystem is formed with the last two lines of system (31) and permits to compute -022 and 
ipi2 with the same procedure. It is interesting to remark that for these cases, the above process always gives 
<Pi2(z) = <f2i(z) = 0. 

3.6. Multilayered angle-ply othotropic plates 

The above procedure is applied. It gives non null tpi2 and <p2i functions. The choice for 7 13 15 etc. is not 
as obvious as for the cross-ply case. The functions that are built with this process seem to have good shapes, but 
it is difficult to verify them with finite element simulations because boundary conditions are difficult to choose. 
All the previous finite element computations were cylindrical bending of the plate in the x or y directions. For 
off-axis plates, it is always possible to enforce cylindrical bending but coupled effects makes the results difficult 
to interpret. Further studies are needed to know what is the better way to determine warping functions in this 
case. 

4. Results 

4-.1. The homogeneous -single layer- case 

Here are presented the warping functions for the homogeneous case. We shall distinguish the isotropic and 
the orthotropic case for two reasons. The most evident reason is that, for an orthotropic material which is 
considered off-axis, four non-null warping functions are expected, instead of only two (different) for the on-axis 
configuration, and one (two identical) for the isotropic case. But there is another reason, also very interesting. 
For this material, it is possible to give to the shear transverse modulus G13 and G23 arbitrary (small) values. 
We are going to see that this has an influence on the shape of the warping functions. 




CV?1 + sV°2 
SV?1 + cV°2 
-CSlfii! + SClf22 
-CSlfi° n + SC(fi22 



(43) 
(44) 
(45) 
(46) 




^21 
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4.. 1.1. The isotropic single layer case 

An isotropic material with v — 0.3 is considered here. As we shall see in figure 2, the thickness to length ratio 
h/L has an influence on the warping functions shapes. However Reddy's formula is a very good approximation 
until this ratio reaches the value of 4.1CP 1 which is uncommon for plates studies because three dimensional 
effects generally occurs with such structures. It is the reason why Reddy's model has been extensively used and 
is found to give good results for many practical applications. But we shall see below that the sensitivity of the 
h/L ratio is higher when orthotropic material are considered. 




Figure 2: Warping functions for the case of a single layer isotropic laminate, with v = 0.3 and h/L varying from 2.10" 1 to 8.10" 1 
compared to Reddy's formula 1 — A(z/h) 2 . The two first curves are superimposed. For clarity, the 2D FEM simulation has been 
plotted only for the two higher values of h/L. 



4-.1.2. The orthotropic single layer case 

Ratio h/L varies: The same study is performed for an orthotropic on-axis material for which the ratio 
G\s/Ex has been set to 40, which is a common value for an unidirectional composite ply. It shall be seen in 
figure 3 that, as mentioned above, the h/L ratio has more influence than for an isotropic material. 

Ratio G\ 3 /E\ varies: Now, an orthotropic on-axis material is studied for various values of ratio G13/-E1 
(see figure 4). Other properties have been fixed to: £2 — £3 = E\, V23 = v\ 3 = i/ 12 = 0.3, G23 = G13, 
G12 — £1/2(1 + 1/12) and h/L = 0.1. Figure 4 shows the warping functions, compared to 3D finite element 
computations. 

Off-axis: In this case, the choice of parameters 7° 3 , 723, 7^3 1 , etc. is determinant. Taking values presented 
in the end of the section 3.2 gives warping functions presented in figure 5. It is difficult to propose a finite 
element validation because a simple bending configuration have not yet been found. 

4-. 2. The multilayered case 

4-2.1. The isotropic multi layer case 

In this section, a n-layer material with isotropic plies is studied. Young's modulus have been fixed accord- 
ing to the rule £ 0( jd layer — l/25£ cven i a yer, except for the case n = 3 for which the reciprocal £ e ven layer = 
l/25£ dd layer configuration is also given. Poisson's ratio is 0.3. The thickness to length ratio has been set to 
h/L = 0.1. We shall see in figure 6 the warping functions and their derivatives. 

4-2.2. The anisotropic multi layer case 

In this section, warping functions of some classical laminates are computed. Considered mechanical prop- 
erties of the ply are: E 1 = 400000, £ 2 = £3 = 10000, is 2 3 = viz = v x % = 0.25, G 23 = 60 00, G 13 = G 12 = 5000 
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Figure 3: Warping functions for the case of a single layer orthotropic laminate with G13/E1 = 40 and h/L varying from 10 -3 to 
4.10 - 1 compared to Reddy's formula 1 — 4(z/h) 2 . The three first curves are superimposed, and the fourth one only differs from 
Reddy's formula of about 0.02%. For clarity, the 3D FEM simulation has been plotted only for the four higher values of h/L. 

The [0/90] laminate: Unfortunately, at this time, the above procedure does not give correct warping 
functions. The problem could be due to coupling effects (this laminate has a membrane-bending coupling) but 
it is also the case for the isotropic two-layered laminate of section 4.2.1 for which warping functions were found 
and agreed very well. The problem could also be due to the discontinuity of the ratio between transverse shear 
and longitudinal modulus at the neutral plane. Further investigations need to be done for this configuration. 

The [0/90] s laminate: In this case, the cylindrical bending gives a decoupled strain state, and finite element 
results are in perfect agreement with warping functions computed from the present procedure, which can be 
seen in figure 7. 

The [30/ — 30] laminate: It this case, the procedure gives warping functions which are plausible but no 
finite element proof can been given at this time for the same reasons than those of the [30] one layer laminate 
of section 4.1.2. They are presented in figure 8. 

The [30/ — 30] s laminate: The same remark as above applies. Results are presented in figure 9. 

5. Conclusion 

Transverse shear warping functions for an equivalent single layer plate model have been built from a vari- 
ational approach. Based on a plate model which has been previously used by other authors, the part of the 
strain energy which involves the shear phenomenon is expressed in function of the warping functions and their 
derivatives. Then the variational calculus leads to a differential system of equations which warping functions 
must verify. 

However, solving this system requires the choice of values for the (global) shear strains and their derivatives. 
This is a crucial point because the chosen values may change considerably the shapes of the warping functions. 
For cross-ply laminates, a particular choice, whose justification is given, lead to excellent results in terms of 
the warping functions. For angle ply laminates, the above choice cannot be justified, but at this time, no other 
choice has been found relevant. 

For single layer isotropic and orthotropic plates, an analytical expression of the warping functions is given. 
They involve hyperbolic trigonometric functions. They differ from the z — 4/3z 3 Reddy's formula which has 
been found to be a limit of present warping functions for isotropic and moderately thick plates. When the h/L 
and/or the G13/E1 ratios significantly differ from those of isotropic and moderately thick plates, a difference 
between present warping functions and Reddy's formula can be observed, even for the isotropic single layer 
plate. Finite element simulations agree perfectly with the present warping functions in these cases. 
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Figure 4: Warping functions for the on-axis orthotropic single layer case with G13 varying from Ei/10 to Ei/10 4 . 

For multilayer cross-ply configurations, the warping functions are expressed from the differential system in 
the form of layer-wise analytical solutions involving matrix exponentials and integration constants which must 
be determined. Warping functions and their derivatives must verify conditions at the bottom, the middle, and 
the top of the plate and at each interface. These conditions form a last linear system whose unknowns are the 
expected constants. Present warping functions have been compared to results of 3D finite clement simulations. 
They are in excellent agreement. 

For angle-ply laminates, the above process gives warping functions that seem to have relevant shapes, even 
if the choice for global shear values cannot be justified, as it was told before. In addition, no finite element 
comparison has been presented because it is difficult to propose boundary conditions and prescribed load that 
permit to isolate the shear phenomenon. Further work must be done to clarify this point. 
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